Statistical ensembles of virialized halo matter density profiles 
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We define and study statistical ensembles of matter density profiles describing spherically symmet- 
ric, virialized dark matter haloes of finite extent with a given mass and total gravitational potential 
energy. We provide an exact solution for the grand canonical partition functional, and show its 
equivalence to that of the microcanonical ensemble. We obtain analytically the mean profiles that 
correspond to an overwhelming majority of micro-states. All such profiles have an infinitely deep 
potential well, with the singular isothermal sphere arising in the infinite temperature limit. Systems 
with virial radius larger than gravitational radius exhibit a localization of a finite fraction of the 
energy in the very center. The universal logarithmic inner slope of unity of the NFW haloes is 
predicted at any mass and energy if an upper bound is set to the maximal depth of the potential 
well. In this case, the statistically favored mean profiles compare well to the NFW profiles. For very 
massive haloes the agreement becomes exact. 



I. PRELIMINARIES 

Numerical simulations of structure formation in the 
current successful cosmological paradigm of cold dark 
matter exhibit long-lived virialized structures, dark 
matter haloes. Their spherically averaged density 
profiles as function of radius are well described by a 
single parameter family of profiles (NFW) [H-Q- Sec 
Q for a recent review and more references. Moreover, 
these profiles are universal with respect to the initial 
conditions of the simulations, yet, it is still unclear to 
what extent this universality can be explained from 
first principles. In particular, the interplay and relative 
importance of statistics and gravitational physics are 
not fully understood. Motivated by these results and 
unsolved issues, we pose a relatively simple problem: 
given only the mass, total gravitational potential energy 
and spherical symmetry of an ensemble of virialized 
haloes, we seek to find the functional form of the most 
likely density profiles. 

The present work is not the first attempt at the 
description of self-gravitating astrophysical systems with 
statistical arguments. Such ideas started with violent re- 
laxation and the classic paper (LB), where the phase 
space distribution f(x, v) maximizing Shannon's entropy 
at a given energy is derived [see also H, for a recent 
extension of that approach directed to similar aims]. The 
resulting phase space distribution is that of an isother- 
mal sphere, that does not reproduce properly the profiles 
found in the simulations. Unlike LB, the present paper 
is not concerned with the phase space distribution of a 
single halo, but studies statistical ensembles of haloes. 
Next we discuss our perspective and its connection to LB. 

Imagine an ensemble of approximately spherically 



symmetric, virialized haloes consisting of N particles 
of some mass m, and of a finite extent r vir , the virial 
radius. If all we are interested in are the spatial degrees 
of freedom, the complete statistical description of this 
ensemble is provided by the joint probability density 
Pn(xi, ■ ■ ■ ,Xn) of observing particle 1 at xt, particle 2 
at X2, etc. According to the virial theorem Q, the total 
energy of each halo is half its gravitational potential 
energy W. A tentative description of the halo ensemble 
may be for example through a canonical ensemble 
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The function pjy is the distribution corresponding to 
maximal entropy for the prescribed mean value of the 
gravitational energy. The study of the associated par- 
tition function and typical realizations does not appear 
simple in this discrete formulation. Neither it is obvious 
how to enforce spherical symmetry in (TTJ) . The principal 
purpose of this paper is to solve exactly the statistical 
ensembles analog to ([1]), where the probability density 
Pn(xi,--- , xjv) is replaced by a probability functional 
p[p(x)], and p(x) is a spherically symmetric positive 
function. 

The non-consideration of the velocity degrees of freedom 
and thus of the kinetic energy is far from innocent. 
Together with the presence of a sharp boundary r v - lr , 
it is in fact one of the elements making our statistical 
ensembles well defined and leading to meaningful results. 
It is indeed well known that the thermodynamics of grav- 
ity is intrinsically difficult 0-[H|; and m particular the 
statistical ensembles are often not well defined due to the 
unbounded available phase space. Identifying the energy 
of the halo with its gravitational energy with through the 
virial theorem short-cuts many such difficulties. On the 
other hand, the specific thermodynamics associated with 
our ensembles are generally different from that of gravity. 

The connection to LB can be made apparent with 
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the following argument. The probability pi(x) to ob- 
serve a particle (particle 1, or any other) at x is obtained 
by integrating out the N — 1 other coordinates of pn- 
We enforce spherical symmetry by enforcing the same 
on pi. A realization of the density within the ensemble 
is p(x) = J2i $ D {x — Xi), and the mean density profile is 
directly proportional to p±, 

( P (x)) =N Pl (x). (2) 

Now, let us assume hypothetically that the following 
uncorrelated form of the N particle probability density 
holds : 

N 

p N {x u - ■ ■ ,x N ) = \\_p{xi). (3) 

i=l 

The entropy of pm is now simply N times that of p\. 
Maximizing the entropy with fixed mean gravitational 
potential energy leads now to 

pi(x) oc (p{x)) oce^ m <*^>. (4) 

This is identical to the density profile obtained by LB, 
isothermal spheres with profiles decaying as p(r) oc r~ 2 . 
In this equation, (<$>(x)) is the gravitational potential 
at x, obtained consistently from the mass distribution 
(p(x)) . In the LB approach, (3 is the inverse temperature 
of the velocity distribution, a Maxwell-Boltzmann 
distribution. 

Thus, one interpretation of the approach taken in 
this paper, especially for the grand canonical ensem- 
ble, is the LB entropy maximization approach to the 
A^-particle distribution in a continuous description, 
ignoring the velocity degrees of freedom. Indeed, our 
approach recovers the singular isothermal spheres in a 
particular regime, but also interesting new phenomena 
arise such as localization of energy, or the possibility of 
negative temperatures. 

This paper is organized as follows. In section |TT] 
we solve the statistical ensembles associated to our 
problem. More specifically, we first find the partition 
function of the grand canonical ensemble in section 
III A[ and then discuss in III Bl its equivalence to the 
microcanonical ensemble in the relevant limit. We 
obtain the mean density profiles that can be realized in 
the most micro-states. We then compare these profiles 
to NFW profiles in section IIIII We conclude with a 
discussion in section HVl 

II. SOLVING THE STATISTICAL ENSEMBLES 

For a spherically symmetric halo of extent r v - 1T with 
density profile /o(r), the mass and gravitational energy 
are given by 

M[p] = Air / drr 2 p(r), (5) 
Jo 



and 




where M(< r) is the mass enclosed within r. The micro- 
canonical ensemble, expressing uniform statistical weight 
on the surface of constant mass and energy, is formally 

fl(M, W) cx J Vp 6 D (M - M[p])5 D {W - W[p\) (7) 

where the integration runs over all positive functions p(r) 
in a sense that we will make precise later, and S D is 
the Dirac delta function. Similarly, the grand canonical 
ensemble reads as 

Z[\,(3}<x Jvpe- XM W-P™. (8) 

We present a solution for the latter first, establishing 
their equivalence afterwards. The parameter can 
be interpreted as an inverse temperature I/HbT, and 
the parameter A as —p/ksT, where p is the chemical 
potential. Equivalently, according to Jaynes' [H, [l3[ 
perspective on statistical mechanics, the corresponding 
probability density p[p] oc e ~ XM ~P\ w \ is the probability 
density of maximal entropy at a given mean mass and 
energy. 

Before calculating the partition function, let us in- 
troduce the following function of radius 

y(r) := Air J ds sp(s) (9) 

This function turns out to the the most fundamental for 
our purposes, as it allows the factorization of the parti- 
tion function. Another representation of y(r) is 

d>(r)=-- f dsy{s), (10) 
r Jo 

where $(r) is the gravitational potential at r, thus the 
average of y(r) between and r. From its very definition 
the y(r) function is always a positive, decreasing function 
of radius until it reaches zero at the virial radius. Further, 
the maximal value y(0) that the y-function takes is the 
central potential of the mass distribution. From equation 

m 

*(0) = -Gi/(0). (11) 
Note that both sides of that equation can be infinite. 

The y-function has the following two convenient 
relations: for any spherically symmetric distribution 
holds that the mass and energy are its first and second 
moments: 

M= f dry{r) (12) 
Jo 
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as well as 



and the dimcnsionlcss density profile 



W = ~jJ dry 2 (r). 



(13) 



The relevant, dimensionless parameter in this work is the 
ratio a of the virial radius to the gravitational radius 



where 
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(14) 



(15) 



For future use, we introduce the depth of the central 
potential well as another parameter 



1 < 



1 $(0) 



e$ ' $(r vir ) 



«<°> (if) 



(16) 



One can show from equations (|12l) and (|13l) that any mass 
distribution must have 



1 1 

- < a < . 

2 ~ ~ 2e$ 



(17) 



The singular isothermal sphere p oc r 2 has a = 1 , e$ = 
and plays a special role in this work. 



A. The grand canonical ensemble 

In terms of the y- function, the partition function of the 
grand canonical ensemble is given by 

Z[X,0] = j Vy{r)e- X ^" dry{r] ~^ dry ^ r \ (18) 



with boundary conditions 

2/(rw) = 0, y{r) > 0, and y'(r) > 0. (19) 

The Jacobian of the transformation from p to y is of no 
relevance, since the relation is linear. The parameters 
and A are related to the mean energy and mass by the 
two non-linear equations 



\W\ = — 



d\nZ 



and 



M = — 



dlnZ 
dX ' 



(20) 



(21) 



From here on we work with dimensionless functions and 
potentials A, 0. We define the dimensionless y-function 



zO) = (^) y(r), x 



r 

: 

T* v ir 



(22) 



u(x) 



z'{x) 



^ [ # ) P(r)- 



(23) 



We approach the formal expression (|18[) through dis- 
cretization, replacing the continuous function by its 
sample on a regular grid of N points, replacing integrals 
by finite sums, and derivatives by finite differences. We 
then obtain a solution that has a well defined limit for 
N -> oo. 

We define Zi := z(xi),Xi := i/(N+l), i = 1, - ■ ■ ,N, 

and set further J Q dx z(x) w jjJ2k z k an d 

Jq dx z 2 {x) w ^Sfc 2 ?- After renormalization of 
the potentials, the partition function is 

Z[X, 0, N] = J °° d N z Uti CX P (-^ (24) 



zi> z 2 > 



> 



■N- 



(25) 



Using the symmetry of the integrand, this is simply 
Z[\,P,N] = -^Q dz cxp(-Az-/3z 2 )^) 



(26) 



The defining equations for A and j3 are seen to become 
dhxZi 1 



and 



1 



2a = - 



= _ / dz z e~ Xz - f3z 



y 2„-\z-Pz 2 



dz 



00 Zi J a 



(27) 



(28) 



Interestingly, A and are now set by the simple matching 
of the first two moments of the one-dimensional distribu- 
tion 



p(z) = — exp (-Xz - 0z 2 ) , 
6\ 



(29) 



a Gaussian distribution with range truncated to 
< z < oo. 

This system of equation has a solution for 1/2 < a < 1, 
that can easily be obtained with standard numerical 
methods. The solution is shown on figure [T] The 
lower bound simply expresses the fact that there are no 
spherically symmetric distribution with a < 1/2. The 
upper bound a = 1, corresponds to infinite temperature 
= 0, and shows a transition to a regime that the grand 
canonical ensemble fails to describe properly. We will 
solve this regime later in section Hi B[ 

The probability density we encounter in (1241) is 
very known from the field of order statistics |l4l. [l5j. 
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involves at any finite N the finite differences 
Ui oc N(z.i — Zi + i)/xi. Fluctuations around the 
mean profile ([55]) do not decay in the large N limit. 
Rather, it can be shown that the density profile at dif- 
ferent arguments behave like a collection of independent, 
exponentially distributed variables centered on the mean 
profile ([55)1 . This is clear for instance for the case of 
p = where the ensemble oc e~ AM does not introduce 
correlations between the density at different arguments. 
This is illustrated on figure [21 showing in the upper 
panel a realization of the function z(x) and in the lower 
panel that of the the density profile, for a = 0.7 and 
N = 1000. The black lines are the exact mean values in 
the large N limit, equations ([52)) and (J55J. 



FIG. 1. The dimensionless parameters A and /3 of the grand 
canonical ensemble as function of r v ir/V s = a. These parame- 
ters are set by equations (|27|l and (|28[) . These equations have 
no solution for a > 1. This regime is better understood with 
the microcanonical ensemble. 



It is the probability density describing the ordering 
in decreasing order of N values drawn independently 
from the same distribution p{z). The properties of the 
micro-states z(x) generated by the ensemble that we 
expose in the following can therefore be obtained from 
well known results with no difficulty. Also, micro-states 
can be generated easily following that procedure. 
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In particular, the ensemble makes for large N a 

sharp prediction on the value of z(xi), fluctuations 

decaying to zero for large N. For < z < oo, let F(z) 
be 



F{z) = ~kl dye ~ XV ~ Py2 
and Q(x) be its inverse function 



< F < 1 



Q(x) = F- 1 (x). 



(30) 



(31) 



In statistical jargon, F is the cumulative distribution 
function of p(z) and Q the quantile function. We have 
that the z-function singled out by the ensemble is 



%i = Q (1 - Xi) . 



(32) 



In the continuous limit, z(x) = Q(l — x). The corre- 
sponding density profile is given by ([23)1 . with the result 



u(x) = cxp (lnZi[A, P) + Xz(x) + /3z 2 (x) -Inx). (33) 

The mass and energy of the micro-states also become 
sharply defined for large N, since they are obtained inte- 
grating z(x). This fact already points to the equivalence 
of the grand canonical and microcanonical ensembles 
discussed later in more detail. 

On the other hand, the density profile u(x) = z'{x)/x 



■ 0000 



1000 r 




FIG. 2. A realization of a micro-state of the grand canonical 
ensemble pi]) (dotted lines), for a = 0.7 and N = 1000, to- 
gether with their exact mean value (solid lines). The upper 
panel shows the dimensionless y-function. The lower panel 
displays the dimensionless density profile u(x). For large 
N, the ensemble always makes a sharp prediction on the y- 
function and other integrated quantities, but realizations of 
the density profiles themselves fluctuate around the mean pro- 
file. 



For j3 — 0, corresponding to the upper bound 
a = 1, we have further through direct inspection 
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Q(x) = — ln(l — x), so that 

z(x) = -In (a;) , /3 = (34) 

and we recover thus the profile p ot r~ 2 of the singular 
isothermal sphere. At finite temperature, an exact if not 
necessarily very insightful form for z(x) is 

z{x) = —= [erfrT 1 (x erfc (7)) - 7] , /3 > (35) 
VP 

with 7 = X/2^/P, and where erfc(x) is the complemen- 
tary error function. Finally, the lower limit a = 1/2 
corresponds to p(z) ps S d (x — 1), leading to the degener- 
ate profile where all the mass is concentrated at the virial 
radius. The outer slope is found to be 

d In p _ 2 

— : — = -1 - Vtt7 e 7 erfc(7), (36 

dmr r=rv . r 

and transitions smoothly from +00 at a = 1/2 to —2 at 
a = 1. 

We now turn to the more careful study of the regime 
a > 1. 

B. The regime r v i r /r g > 1 and the microcanonical 
ensemble 

In the discretization introduced above, the micro- 
canonical ensemble takes the following form, 

n N = J d N z 6 D (1 - i E k z k ) 5 D (2a - i £ fc $ 

00 > z\ > z 2 > ■ ■ ■ > z N > 0. (37) 

By symmetry, it reduces to 

n - = m r ^ ( x - * ? s ° ( 2a - * ? **) 

. (38) 

The mass (in our units unity) defines a plane in N- 
dimcnsional cartesian space, while the energy 2a defines 
a sphere. The microcanonical ensemble is remarkably 
reduced to the (rather formidable) geometric problem of 
the intersection of a plane and a sphere in N dimensions, 
together with the constraints z, > 0. A realization of z 
is obtained by drawing at random a vector on this set, 
and then ordering its coordinates. 

Unlike for the grand canonical ensemble, the calcu- 
lation of the y-function at fixed N does not appear 
tractable with simple means. Nonetheless, the uniform 
distribution on this set (without ordering) has already 
been studied in the limit of large N [lj| , and we can use 
the main results exposed there. Notably, the equivalence 
to the grand canonical ensemble for 1/2 < a < 1 rigor- 
ously follows directly from Theorem 1.1 in that work jl7| . 



Theorem 1.2 of [l6| deals with the case a > 1, 
demonstrating localization of energy. It is shown that 
the largest component of such a random vector z (in 
our notation z(l/N)) grows in such a way with N that 
the part of energy -^z 2 (l/iV) that it carries is finite 
with probability 1. It is not a localization of mass in 
the sense that the fraction of mass j^z(l/N) still tends 
to zero. The amount of energy localized in the center 
is precisely the surplus of energy with respect to the 
singular isothermal sphere a = 1. The second largest 
coordinate does not show this localization property. 
The density profile that we recover is in fact exactly 
that of the singular isothermal sphere with a stronger 
singularity in the very center x = accounting for the 
additional energy. Note that the situation is reminiscent 
to Bose condensation. 

The above results can also be seen from the per- 
spective of a modified grand canonical ensemble. The 
grand canonical ensemble fails for a > 1 because a 
positive inverse temperature j3 is required, for the 
system can reach arbitrarily high energies \W\. With 
j3 < the probability density p(z) cannot be normalized 
anymore. As in other physical systems we can allow 
negative temperatures by setting an upper bound to 
the energy that the system can have access to. This 
is simply done by introducing a maximal value for the 
z-function. Since z(0) is the potential well l/e$ in (fi"6|) . 
this is equivalent to force the central potential of the 
mass distribution to take a finite value. 

Our results of the previous section, equations (|24|) 
to (f3U| all hold unchanged provided the range of z in 
p(z) is now < z < l/e$. In particular, to each a 
within the allowed range 1/2 < a < l/2e$ there is 
a corresponding A and /3. Since z(x) is given by the 
quantile function Q(l — x) of p{z) we have 



z(0) = -. (39) 



The profile singled out by the ensemble always saturates 
the bound. 

To very small e$ and a > 1 correspond now high 
negative temperatures /3 <~ and A w 1. This 
implies that p(z) is over a long range the exponential 
distribution, with relevant deviations only for z w l/e$ 
accounting for the match of the second moment of p(z). 
Thus we recover as we just argued the exact singular 
isothermal spheres with the corresponding additional 
singularity in the center. As we show next, the limit 
introduced for the central potential has the additional 
benefit of recovering the universal NFW profiles seen in 
simulations. 
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III. NEGATIVE TEMPERATURE PROFILES 
AND NFW HALOES 

The one-parameter family of NFW profiles is given by 



1 



unfw{x) oc — 



1 



x(l + ex)' 



(40) 



where c is the concentration parameter of the halo. More 
massive haloes have smaller concentration, and the range 
of haloes probed in numerical iV-body simulations corre- 
spond approximately to the range 3 ~ c ~ 15. The 
corresponding y-function is given by 



znfw{x) oc 



1 



1 



1 



From this relation follows 



e*,NFW = -J [(1 + c)ln(l 



as well as 



onfw 



J_ ( 1 — 2e$ ; NFW 
2c \ £$.nfw 



(41) 



(42) 



(43) 



It is clear that the ensembles at positive tempera- 
ture (3 > are unable to accommodate for these profiles. 
Our study of these ensembles in the previous section 
showed that the overwhelming majority of profiles at 
given mass and gravitational potential energy have their 
y-function given by (|35l) for a < 1, or basically by the 
singular isothermal sphere (|3"4")l for a > 1. None of these 
profiles can provide a satisfactory match. The outer 
slope at x = 1 is not steep enough, and the singularity 
in the middle is too strong. The NFW density profile 
is singular in the center, but as can be seen from (|4"Tj) 
all NFW haloes have a finite central potential z(0), 
cquivalcntly e$ > 0. The profiles singled out by the 
ensembles all have z(0) = oo, e$ = 0. 



The dashed line in Figure [3] shows the the loca- 
tion of the NFW haloes with 3 < c < 15 in the 
(e$ — a) plane, together with several lines of constant 
(3 = 1,0.5,0,-0.5,-1 (solid lines, from bottom to top). 
The dotted lines a = 1/2 and a = l/2e$ corresponds to 
p = ±oo, or zero temperature, where there is only one 
micro-state for the given parameters. The line (3 = 
converging to a = 1 for small e$ marks the onset of 
negative temperatures, where the number of available 
micro-states starts to decrease with increasing a. The 
parameters of NFW haloes are all assigned negative 
temperatures. Interestingly the NFW haloes follow 
rather closely the line of constant (3 « —0.17, the dotted 
line. 




FIG. 3. The dashed line show the location of the spherically 
symmetric NFW haloes (3 < c < 15) in the plane set by 
a = r v ir/r s and e$ = $(r v i r )/"I > (0), with more massive haloes 
having larger e# . There is no spherically symmetric mass con- 
figuration outside the region delimited by the two dotted lines 
a = 1/2 and a = l/2e*. The solid lines illustrate the behav- 
ior of the line of constant (3 of the grand canonical ensemble. 
The NFW haloes are seen to lie close to the line of constant 
(3 = -0.17, the dotted line. 



If some process is preventing the growth of the central 
potential, at least over the relevant time scales, we 
can still ask whether the ensembles with the additional 
constraint of a finite potential well introduced in section 
III Bl can reproduce these profiles more successfully. 

There is a strong phcnomcnological argument sup- 
porting this view. As noted previously, the y-function of 
these ensembles are smooth functions that now tend to 
the finite value l/e$ at x = 0. Therefore, the behavior 
of the density profile z'(x)/x in the center predicted by 
these ensembles is precisely the power law of the NFW 
profile 



z(x) oc — , x — )■ 

x 



(44) 



at any mass and energy. 



As discussed in III Bl the only modification to the 
case of unconstrained e$ is the restriction of the range 
of p(z) to < z < l/e$, and the corresponding change 
in its quantile function Q. While an exact but lengthy 
and not very insightful expression for z(x) = Q(l — x) 
can still be written in terms of the error function for any 
value of A and /3, we prefer to give a characterization 
through the differential equations 



z'(x) = exp (in Z\ [A, 0\ + Xz(x) + fiz 2 {x)) = u(x)x, 



and 



z"(x) = {z'{x)) (A + 2(3 z{x)) 



The boundary conditions are 

z(l) = 0, z(0) 



c<i> 



(45) 



(46) 



(47) 
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Figure [4] compares the y-function (upper panel, nor- 




0.0 0.2 0.4 0.6 0.8 1.0 

r / r vir 




0.001 0.010 0.100 1.000 

r / r vir 



FIG. 4. Upper panel. The solid lines show the y-functions of 
the NFW profiles, for c — 3,5, 10, 15, from top to bottom. The 
dashed line the y-function singled out by the grand canonical 
ensemble at the same parameters a and e*. Lower panel : 
the same for the density profile, on a logarithmic scale with 
arbitrary normalization. See text for more details. 



malized to unity at x = 0) and the density profiles 
(lower panel, with consistent but otherwise arbitrary 
normalization) of the NFW profiles and the one resulting 
from the statistical ensemble for the same a and e$. 
The agreement is very good both qualitatively and 
quantitatively for massive haloes. At high concentra- 
tion, too much mass is assigned to the transition region, 
and the outer slope is not well reproduced anymore. 
Interestingly, if e$ is progressively pushed to zero from 
this point, the transition region is pushed back toi^O, 
the mass assigned to the transition region becomes 
the additional singularity in the center discussed in 
III Bl and the outer slope becomes uniformly —2, the 
isothermal sphere. The apparition in our ensembles of a 
characteristic scale radius, roughly x = 1/c for the NFW 
profiles, is a direct consequence of the finite value forced 
upon the central potential. 



IV. SUMMARY AND DISCUSSION 

In this paper we defined statistical ensembles of 
virializcd spherical halo density profiles at a given mass 
and gravitational potential energy. More specifically, 
we first discrctized the relevant radial functions of 
radius on a regular grid and then obtained their exact 
thermodynamics in the limit of vanishing cell size. 
We obtained the mean density profile predicted by 
theses ensembles, and showed that the ensembles make 
in fact a very sharp prediction on the y-function of 
the most likely profiles under these constraints. The 
profiles we found in that way all have a infinitely deep 
potential in the center, and we recovered the well 
known singular isothermal sphere in the case where the 
virial radius is equal to the gravitational radius. Sys- 
tems with larger ratio show the interesting phenomena 
of localization of the additional energy in the very center. 

It is plausible, and has already been suggested [H, fl9j . 
that the growth of the central potential is prevented, 
or is very slow, for cosmological haloes. Implementing 
this constraint in the ensembles leads to the prediction 
of a universal inner logarithmic slope of —1, as for the 
NFW profiles seen in cosmological simulations. At the 
corresponding values of mass, energy and depth of the 
potential well, we showed that the NFW profile is very 
close to the mean value of the statistically preferred 
profiles, especially the massive haloes. Due to this ad- 
ditional constraint, these systems are assigned negative 
temperatures. The number of micro-states available is 
actually decreasing if additional energy \W\ is given to 
the system. This suggests that the NFW profiles might 
not be the end stage of gravitational evolution, rather, a 
long lived quasi equilibrium state. 

It should be stressed that the statistical arguments 
presented here do not explain from first principles the 
emergence of these profiles : an explanation for the 
location of the NFW haloes in the relevant parameter 
space shown in figure [3] is left for future research. 
This location might be related to the type of initial 
conditions encountered in cosmological simulations, 
or to any physical process that our approach cannot 
accommodate. Nevertheless, the appearance of the inner 
slope of unity, and that of a characteristic radius is a 
common feature of all the profiles singled out by the 
ensembles at any mass and energy, as long as the central 
potential is forced to take a finite value. These results 
possibly shed new light on the respective roles of physics 
and statistics in the emergence and stability of these 
profiles. 
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